Scale-renormalized matrix-product states for correlated quantum systems 
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A generalization of matrix product states (MPS) is introduced which is suitable for describing 
interacting quantum systems in two and three dimensions. These scale-renormalized matrix-product 
states (SR-MPS) are based on a course-graining of the lattice in which the blocks at each level 
are associated with matrix products that are further transformed (scale renormalized) with other 
matrices before they are assembled to form blocks at the next level. Using variational Monte Carlo 
simulations of the two-dimensional transverse-field Ising model as a test, it is shown that the SR- 
MPS converge much more rapidly with the matrix size than a standard MPS. It is also shown that 
the use of lattice-symmetries speeds up the convergence very significantly. 
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In a variational study of an interacting quantum sys- 
tem, a wave function with a number of adjustable pa- 
rameters pi , . . . , Pra is optimized by minimizing its energy 
E = {^\H\'^) with respect to a haniiltonian H. Ideally, 
one would like to consider a functional form which allows 
for a systematic way of improving the calculation, by in- 
creasing the number of parameters m in such a way that 
^' is guaranteed to approach the true ground state of H in 
the limit m — > oo. A trivial way, in principle, is to expand 
l^*) in a complete set of states; l^") = X]nC„|n), whence 
the parameters to be optimized are the wave function co- 
efficients c„ themselves. However, in practice, the Hilbert 
space is too large (2''^ states in the simplest case of TV 
spins with S = 1/2) to include all states, and in general 
there is no obvious way to order the states so that their 
contributions to the ground state decrease as a function 
of n. To achieve this, one can attempt to optimize the 
basis in some way, with the goal of obtaining a hierarchy 
of basis states which systematically and rapidly improve 
the result as they are included in the calculation. This is 
the basic idea of renormalization group (RG) methods, 
which after decades of attempts, following Wilson's pio- 
neering solution of Kondo impurity problem pj , led to a 
break-through in the form of White's density matrix RG 
(DMRG) method [3,[3 for one-dimensional systems. 

For systems in higher dimensions, there has been re- 
cent progress in generalizing the DMRG approach, which 
is closely related to matrix-product states (MPS) [4,], us- 
ing tensor-network states |5|, e.g., the projected entan- 
gled pair states (PEPS) [6| and related MPS-like string 
states 0, as well as schemes based on entanglement 
renormalization Q. However, there are still consider- 
able challenges related to the convergence properties and 
computational complexity of these methods. In this Let- 
ter, an alternative class of generic correlated states — 
scale-renormalized matrix-product states (SR-MPS) — is 
introduced. These states combine concepts of coarse 
graining, renormalization, and MPS into a framework 
for systematically refined variational calculations. The 



scheme is tested on the two-dimensional transverse-field 
Ising model, using a recently developed variational Monte 
Carlo method to optimize the SR-MPS. 

In the DMRG method, the basis is re-optimized and 
truncated at some number D of states as more sites are 
added to the lattice . Originally this was not viewed as 
a variational method, but it was soon recognized that the 
DMRG in effect produces the best variational MPS with 
D X D matrices [J| . For a system of N spins represented 
by Pauli operators cji , and working in the basis where all 
af are diagonal, af — ±1, MPS for a periodic chain are 
of the form (using the notation [cr] for [af,..., cr^] ) 



where the wave-function coefficient is 

W{[a]) = Ti'{A{at)A{a^,)---A{<j^^)}, 



(1) 



(2) 



and A{±1) are two D x D matrices. For systems with 
open boundaries, for which DMRG and standard MPS 
techniques are best suited in practice, the matrices are 
site dependent, and in stead of taking a trace the edge 
matrices are vectors. The DMRG method does not oper- 
ate with MPS explicitly, but recently schemes have been 
devised for working directly with the MPS without in- 
voking the DMRG procedures. This formally reduces the 
scaling of the computational effort from to for pe- 
riodic systems 10]. Using Monte Carlo sampling of the 
spins states, instead of evaluating their traces exactly, 
the scaling can be further reduced to [tI, 0] . 

In higher dimensions, the DMRG method typically is 
implemented by regarding the system as a chain folded 
up to form the lattice of interest 1]|. In this effective 



one-dimensional system there are long-range interactions. 
Further studies of the MPS formalism also showed why 
the DMRG method performs poorly in this case. The 
exponential scaling in the number of states that has to 
be kept [l^, or, equivalently, the matrix dimension D of 
the MPS, is a consequence of the inability of the matrix 
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FIG. 1: (Color online) SR-MPS hierarchy for a square lattice. 
The circles represent the original matrices B'^^y — A{ax,y)- 
Products of four of these are indicated by connecting lines. 
A square enclosing such a unit represents scale renormaliza- 
tion with matrices M2,M^, which results in B^ y. Succes- 
sive levels of connected squares enclosed by larger squares 
represent products of four matrices B" y followed by scale- 
renormalizations with M£,M^, resulting in B"^^ . 

products to account for entanglement between neighbor- 
ing sites when the corresponding matrices are far apart 
[13|. To circumvent this problem, tensor-network states 
have been proposed as natural and effective generaliza- 
tions of the MPS/DMRG to higher dimensions [1,11]. 

Here SR-MPS is proposed as an alternative general- 
ization of MPS for higher-dimensional systems. For a 
periodic system of S* = 1/2 spins <Tx,y, each lattice site 
{x,y), x,y = I, . . . , L, is associated with a D x D matrix 
A{a^ y) as in the MPS. However, instead of just arrang- 
ing these matrices according to a string on the lattice, 
the system is first subdivided into blocks, which are as- 
sociated with matrix products. These matrix products 
are then scale renormalized by transforming them with 
some other matrices, before they are multiplied by simi- 
lar block-matrices to represent a larger cell. For a square 
lattice, the resulting hierarchy of matrix products (in the 
simplest case based on blocks with four sub-blocks at 
each level) is illustrated in Fig. [T] The block-matrices at 
levels n and n -I- 1 are related according to 



Note that under the trace of the final assembly of 
products, i?i 1, an equivalent way of defining the block 
matrices ([3]) is with a single scale-renormalization ma- 
trix; BJ+i = BJ,^BJ+i,^BJ+i,,,+iS?.,,+iM". However, 
Eq. ([3]) allows for the possibility of increasing the matrix 
size with the level n, using rectangular matrices M£ and 
Af]J of size Dn+i x and £>„ x Dn+i, respectively. 
This may be useful if the scale-renormalization is fur- 
ther refined by making M£ ^ dependent on the physical 
state of the blocks they transform, using, e.g., a block 
spin Tix,y = 0, ±1 [for six different matrices M£ j^{'>lx,y)]- 
Defining an appropriate block-variable for a given model 
is not always easy, however. In the Ising model con- 
sidered here the Kadanov block-spin lj| can be used, 
but in this first study only the state-independent scale- 
renormalization ([3]) will be applied and the matrix size 
will be kept constant; Z?„ — D. Keeping the symmetric 
form with left and right scale- renormalizations, instead 
of just a single M", seems to help in the optimization, in 
spite of the larger number parameters. 

In the simplest version of the SR-MPS the wave func- 
tion coefficient is the trace of _B' = B[^. One can also use 
a sum of matrix products taken over symmetry transfor- 
mations of the spin configuration. Spin-inversion sym- 
metry can be used if the hamiltonian has it. Then 



W{[a])=Tr{B\[a])±B\-[a])}, 



(4) 



where — [a] denotes the configuration with all ct^ — > — cr^. 
It is also useful to incorporate lattice symmetries. De- 
noting a transformation (including the identity) of [a] 
(translation, rotation, or reflection) by rr[cr], the wave 
function coefficient is, considering for simplicity a fully 
symmetric wave function with zero momentum, 

Wi[a]) = J2Tr{B^{Tr[a]) + B^i-TrW])}. (5) 

r 

Here r — 1, . . . , 8A^ if all symmetries of the square lattice 
are used. It will be shown below that the use of symme- 
tries improves the D convergence very significantly. 

To test the SR-MPS scheme, it will be applied next to 
the Ising model in a transverse field; 



(6) 



with the lowest level corresponding to the original spin 
dependent matrices, B^y = A{ax^y), and M2,M^ ac- 
complishing the scale renormalization. At level n the 
block coordinates take the values fc2", k ~ l,...,L/2". 
Thus the lattice size should be a power of 2; L = 2'. 

The purpose of the scale renormalization is to compen- 
sate, as much as possible, for the non-equivalent ways in 
which the four blocks at a given level n are treated in 
the associated product of four matrices. It will be shown 
that the effect indeed is to make the four members of a 
block more uniform in their correlations with each other 
and the rest of the system. 



L + l 



with periodic boundaries {(jf^^i y = erf ^ and 
(Jx.i)- Computations are expected to be the most chal- 
lenging at quantum-critical points; h in the vicinity of 
he ~ 3.044 |l5| will be the main focus here. 

To optimize the wave function, here using general 
(non-symmetric) real matrices ^(±1) and M£,M]J, n — 
1, . . . , Z, the variational Monte Carlo method discussed 
in Ref. [1| is used. The energy derivatives are calculated, 
and based on their signs the matrix elements are updated 
by a random amount, e.g., atj Uij —sign{dE /daij)Srij . 
Here r^j G [0, 1) is a random number and the maximum 
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FIG. 2: (Color online) Nearest-neighbor spin correlations on a 
4x4 lattice at /i = 3 calculated using (a) an MPS with D — 2, 
(b) an SR-MPS with D = 2, and (c) an SR-MPS with D = 8. 
The thinnest and thickest bars correspond, respectively, to a 
— 11% and +11% deviation from the average. The bars at 
the right and upper edges represent the correlations across 
the boundaries (periodic boundary conditions are used). The 
energies -E/N of the states in (a),(b),(c) are 3.1746, 3.1772, 
3.2108. The exact energy for L = 4 is -3.2155081. 

step S is gradually reduced. If this reduction is suffi- 
ciently slow, the converged matrices will correspond to an 
energy minimum. Calculations in one-dimension Q have 
shown that the global energy minimum, which evolves to 
the true ground-state energy with increasing D, can be 
reached with this method at least up to Z) « 50 

If no lattice symmetries are used, the Metropolis prob- 
ability of flipping a spin can be evaluated with cx op- 
erations using the sequential flip scheme of Ref. [9] . How- 
ever, with symmetries incorporated according to Eq. ([5]), 
the spins cannot be sequentially visited in all trans- 
formed conflgurations, and therefore a different scheme 
has to be employed. Organizing partial products in tree- 
structures, one for each lattice transformation, the total 
number of operations required for each spin update is 
oc N\n{N)D'^, where the factor N is due to the num- 
ber of different matrix products and In(A^) comes from 
recalculating one branch of a tree. 

First, an L = 4 lattice will be considered. To demon- 
strate some of the effects of scale-renormalization. Fig. [2] 
shows a plot of the spatial variations in the nearest- 
neighbor correlations (cT^,yCr|+i,j,) and (crlya^y^^), ob- 
tained with and without scale-renormalization. Spin- 
inversion symmetry is taken into account but no lattice 
symmetries are used, i.e., the states are sampled accord- 
ing to the symmetric (-I-) Eq. ([4]). With a small D one 
would then expect to see traces of the particular way the 
blocks are constructed. Without scale-renormalization 
(i.e., M£,M^ = /), the scheme reduces to an MPS cal- 
culation with the matrix product taken along the partic- 
ular "coarse-graining string" used here in the SR-MPS. 
In Fig. [5Ja), obtained with D = 2, it can be seen clearly 
that the correlations arc non-uniform; in particular dif- 
ferent sites within the 2x2 blocks are not equally corre- 
lated with their neighbors. With scale-renormalization, 
Fig. [2jb) shows a significantly reduced non-uniformity 
between the blocks. The energy is also improved. With 
£) = 8 in the SR-MPS, shown in Fig. ^^c), the non- 
uniformity is much reduced and the energy is improved 
considerably. In contrast, an MPS with D = 8 (not 
shown) only marginally improves on the D = 2 result. 
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FIG. 3: (Color online) Relative deviation Ae ^ {E ~ Ed/E 
of the energy Ed of SR-MPS and MPS with D x D matrices 
from the exact energy E for a 4 x 4 system at h = 

Fig.[3]illustrates the convergence of MPS and SR-MPS 
with an without lattice symmetries. Without lattice sym- 
metries, the MPS energy converges extremely slowly and 
it is not possible in practice to obtain the ground state. 
By incorporating the lattice symmetries the convergence 
is substantially improved, however. In the case of the SR- 
MPS, an exponential convergence with D is seen even 
without lattice symmetries, and with these symmetries 
included the convergence is very rapid, with an accu- 
racy of 10~^ reached already at D = 5. Clearly, both 
lattice symmetries and scale renormalization have very 
favorable effects, and when using both of them the con- 
vergence properties seem very encouraging. 

Before moving to larger lattices, a further improvement 
of the SR-MPS is noted: One can use different matrices 
for all the states of a block of spins, instead of just the two 
matrices A(cr| y) for individual spins. Using 2x2 blocks, 
there are 16 matrices A{a1^ y), where crl^ y = 0, 15 la- 
bels the states of the spins a^ y, a^+i^y, <J^^y+i, a^+i^y+i- 
This gives additional flexibility to the wave function, and 
thus a faster D convergence can be expected. The larger 
number of matrices to optimize does not seem to pose 
any difficulties in practice, and in addition roughly four 
times less operations are required for a spin update. Here 
calculations with 2 x 2-block matrices A{0, . . . ,15), will 
be compared with the basic ^(±1) scheme. 

Fig. m shows the D convergence of the energy and the 
magnetization of an 8 x 8 lattice close to the quantum- 
critical point; h = 3.044. Results for comparison, E/N = 
-3.23627(2) and = 0.14073(1), were obtained with 
the stochastic series expansion (SSE) method flBi- Again 
it can be seen that the use of lattice symmetries is crucial 
for achieving good convergence. The convergence is also 
significantly better with the block matrices A{0, . . . , 15) 
than with single-spins matrices A(±l). 

As with MPS or PEPS the accuracy of an SR-MPS 
calculation with fixed D is expected to be lowest at a 
quantum-critical point. This is explicitly demonstrated 
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FIG. 4: (Color online) Relative error of the energy and the 
squared magnetization for an L = 8 system ai h = 3.044, 
using SR-MPS with single-spin A{±1) matrices and 2x2 block 
matrices A{0, . . . , 15), with an without lattice symmetries. 

for the SR-MPS description of the transverse-field Ising 
model in Fig. \5[ using both single-spin and 2x2 block 
matrices with D = 4. With the single-spin matrices, the 
error in the squared magnetization in the neighborhood 
of the critical field is « 10%, and with the 2x2 block ma- 
trices it is « 3%. The accuracy of the SR-MPS increases 
rapidly away from the critical point. 

In Ref. [Tj] string-state calculations for a 10 x 10 lat- 
tice were reported. A magnetization curve exhibiting a 
phase transition was obtained, but the accuracy is ac- 
tually rather poor close to the critical point, with devi- 
ations of more than 50% from the exact result and too 
little finite-size rounding. Since no systematic conver- 
gence tests were presented it is difficult to compare the 
performance of SR-MPS and string states directly. 

Although the 8x8 lattice considered here is small in the 
context of QMC simulations of sign-problem-free models, 
the calculations demonstrate that the SR-MPS approach 
is a practically feasible. One important question is of 
course how the D required to obtain a desired accuracy 
grows as N increases. In order for the scaling to be a 
power-law, instead of exponential, it is believed that an 
area law for the entanglement entropy has to be obeyed 
17| . Because of the lattice symmetries incorporated, the 
SR-MPS may satisfy such a law [l^. This, however, 
would still not guarantee a power-law scaling. Numeri- 
cally it is also currently difficult to establish the scaling, 
because a range of system sizes are needed. Calculations 
for L — 16 aX h = he give an energy error A^; < 0.2% 
for 13 = 5 . Thus it is at least clear that one can access 
practically useful lattice sizes. 

The SR-MPS scheme also work for frustrated systems. 




FIG. 5: (Color online) Squared magnetization as a function 
of the external field for an L = 8 system obtained with _D = 4 
SR-MPS, using matrices either for single spins of for blocks of 
2x2 spins at the lowest level. The solid curve shows results 
obtained with the approximation-free SSE method. Relative 
errors are shown in the inset. 

Preliminary calculations for an L ~ 8 square-lattice 
S" = 1/2 Heisenberg model with a ratio J2/J1 = 0.5 of 
the second-nearest to nearest-neighbor interaction show 
a slower convergence with D than in Fig. |31 but it does 
appear feasible to reach the ground state. 

It was shown here that incorporation of lattice sym- 
metries in the wave function is crucial for achieving good 
convergence. This has also been noted for one dimen- 
sional MPS [1^. Although the scaling of the computation 
increases by a factor of N, to oc D^N^ ln(7V) operations 
for updating the whole system, this can be partially al- 
leviated by parallelizing the calculation of the spin-flip 
probability — the 8A^ different traces can be calculated 
completely independently of each other. 
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